    

    ## negative log likelihood function
    nLL = function(par.update, parn = "not", pars, data, eta.coef, 
                   cnt, rep=500, seed){
        
        if(parn %in% names(pars)){
            pars[[parn]][1:length(pars[[parn]])] = par.update
        } 
        
        L0 = cbind(rep(1,nrow(data)),data[,1:pL0])
        l0 = L0 %*% 2^(1:(pL0+1))
        L0 = unique(L0,margin = 1)
        l0.index = match(l0, unique(l0))
        Nb = nrow(L0)
        
        theta.coef = matrix(pars$theta.coef[c(1,3:(pL0+2), 2:(pL0+2))],,2)
        colnames(theta.coef) = c("l0","l1")
        theta = exp(L0 %*% theta.coef)
        
        phi.coef   = matrix(pars$phi.coef[c(1,4:(pL0+3),2,4:(pL0+3),
                                            3:(pL0+3))],,3)
        colnames(phi.coef)   = c("a0","a1","a.") 
        phi   = exp(L0 %*% phi.coef)
        eta   = expit(L0 %*% eta.coef)
        gop   = exp(L0 %*% pars$gop.coef)
        
        ### Step 1: get ratio between two "adjacent" probabilities
        ratio = t(sapply(1:Nb, function(i) GetRatio(theta[i,], phi[i,],eta[i,])))
        dim(ratio) = c(Nb, 5, 2)
        dimnames(ratio) = list(1:Nb, c("ak","lk","aK","lK","l0"),c("0","1"))
        
        ### Step 2: get the maximal ratio against the "baseline p000" 
        #   for each observed unit (only depends on L0)
        k.max = GetMaxKGen(ratio, K, Nb)
        k.max[k.max == Inf] = .Machine$double.xmax
        check = which(k.max < 1e10)
        
        ### Step 3: get the largest p (only depends on L0)
        p.max = rep(1 - 0.1^10, Nb)
        if(length(check)>0)     p.max[check] = PMaxEstRough(ratio[check,,], K, k.max[check], gop[check],rep = rep, seed = seed)
        
        if(parn=="n01"){
            return(sum(p.max > 2*.Machine$double.eps & 
                    p.max < 1-0.1^10 - 2*.Machine$double.eps))
        }
        
        
        ### Step 3: get the ratio against the largest p for each observed unit
        ratio = ratio[l0.index,,]
        p.max = p.max[l0.index]
        k.max = k.max[l0.index]
        A = data[,(pL0+1):(pL0+K+1)];   colnames(A) = 0:K
        L = data[,(pL0+K+2):(pL0+2*K + 2)];   colnames(L) = 0:K
        k = GetKGen(A, L, ratio, k.max, K)
        k[k.max == .Machine$double.xmax] = 0
        
        p.y   = p.max * k
        
        Y     = data[,ncol(data)]
        ll    = sum(cnt[Y==0] * log(1-p.y[Y==0])) + sum(cnt[Y==1] *log(p.y[Y==1]))
        
        return(-ll)
    }
    
    
    # system.time({ for(i in 1:100)  nLL(pars, data, eta.coef)  })
    
    
    